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ABSTRACT 9 


This thesis is concerned with a comparative study of discrete 
time filters using the theories of Wiener-Kolmogorov, Bode-Shannon, 
and Kalman, applied to the filtering of a non-stationary random 
signal in the presence of measurement noise. Programs are developed 
for the simulation of these systems and signals on a digital computer. 
Their filtering properties are compared for a random input signal 
with known steady state characteristics, starting at initial time 
t = с Тһе results show that when the gain of the Kalman filter 
is equal to its steady state value, the Kalman and Wiener-Kolmogorov 
filters perform identically. For large initial errors the Kalman 
filter, with large initial gain, gives the best transient response; 
for small initial errors the Kalman and Wiener-Kolmogorov filters 


are essentially equivalent in their transient responses. 
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CHAPTER I 
INTRODUCTION 
1. OBJECTIVE 


This thesis examines and compares several techniques for the 
filtering of random signals in the presence of measurement noise. 
In diagram form, the problem is as shown in Figure 1; x(t) is the 
random signal to be estimated, v(t) is a random measurement noise, 
$(t) is the estimated or filtered signal, and e(t) is the error in 
estimation. The criteria for the optimum filter is one which 
reduces the mean square error to a minimum. The following techniques 


are studied and compared with a common example: 


(1) Discrete Wiener-Kolmogorov filter 

(2) Discrete Kalman filter 

(3) Bode-Shannon discrete time filter 

For all of the filters it is assumed that the statistical 
Properties of the signal and the noise are known by their correlation 
functions, or their power density spectra, or by equivalent white 
noise excited dynamic models. It is also assumed that the signal 
and the measurement noise are uncorrelated. Wiener-Kolmogorov theory 
assumes that the signal and the noise are stationary; that is, their 
statistical properties do not vary with time. This implies that all 
signals are initiated at t = -o, Kalman and Bode-Shannon filtering, 
on the other hand, &ssume that signals start at some initial time, 
tg» and that as time progresses the stationary probability charac- 


teristics of the signals are established. In addition, Kalman (as 


10 


1] 21013 5 


(3)x 


as tou 


yusawsınseoaw 





1/1 120111 221114 31068 


І 3 


URSI IIS 


+ 939ur4se 
(zy 








10119 


UOTIPUTISO 


well as Bode-Shannon) assumes that the random signal process is 
Markovian, that is, can be thought of as generated by a linear 
dynamic system excited by white node cw During the transient 

stage, from t - to until the statistical properties are established, 
the Kalman and Bode-Shannon filters behave as time-varying devices. 
In the steady state, all of the filters can be expected to behave 
identically since they each establish a minimum mean square error 
estimate. 

Τη order to test and compare the operation of these filters, 
particularly for a non-stationary input, the Wiener-Kolmogorov and 
Kalman filters have been designed for discrete time simulation on 
a digital computer for filtering the same random input signal and 
measurement noise. The output of each filter and the actual random 
input signal have been compared starting at time t = to until steady 
state conditions are established. 

The type of filtering problem studied can be compared with a 
basic tracking problem where a target, with statiscally known 
maneuvering properties, suddenly appears. It is desired A track 
the target and, in particular, to investigate how long it takes the 
output to reach the steady state minimum mean squared error. The 
results of this investigation are presented in this thesis together 
with a summary of pertinent background theory. In addition, the 


discrete time formulations for the Wiener-Kolmogorov theory have 


been developed and programmed for the specific example treated. 


*R. E. Kalman, New Methods and Results in Linear Prediction 
and Estimation Theory, Technical Report No. 61-1, (Baltimore, 
Maryland: Research Institute for Advanced Study, 1961) p. 3. 
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Ihe basic theory for discrete Kalman and Bode-Shannon Lilterimp 
appears in the literature and has been adapted to the example. 
Programs necessary for generating the stochastic signals from 


white noise have also been developed. 
2. HISTORICAL BACKGROUND 


The many present day theories for the filtering of time 
sequences are basically derived from the original work of Wiener 


TEL They obtained 


and Kolmogorov published in 1941 and 1942, 
the direct solution of the filter problem independently and 
simultaneously by the use of the calculus of variations. This 
approach led to the filter's impulse response being characterized 
by the solution of a difficult integral equation known as the 
Wiener-Hopf equation. 

One of the next major developments in the solution of the 


2d in 1950 


filtering problem was accomplished by Bode and Shannon 
with the simplification of the derivation of the Wiener-Kolmogorov 
filter. They used frequency domain analysis and conventional 
circuit theory concepts to interpret mathematical operations in 
physically intuitive terms. Bode and Shannon also applied this 
frequency domain technique to discrete signals with constant 
sampling intervals, 

In most cases, the solution of the Wiener-Hopf integral equation 
is a formidable task. With the advent of the digital computer, an 
alternate approach was sought to avoid this difficulty. In 1960 


1ο] 


Kalman and men attacked the filtering problem from the 


state variable point of view. They evolved a set of differential 


12 


equations, equivalent to the Wiener-Hopf integral equation, whose 
solution yields the optimum filter. These differential equations 
may be synthesized in a sequential fashion and therefore digital 


computer solutions easily obtained. 
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CHAPTER II 
IHE WIENER-KOLMOGOROV FILTER 


There are several methods for deriving the continuous Wiener- 
Kolmogorov filter. Appendix À summarizes the solution which leads 


to the Wiener-Hopf integral equation. 
со 


P p (t) = eo P(t - T)dr, t >2 O (2-15 


=0 

H(T) is the desired filter impulse response, 9, E) is' the auto- 
correlation function of the signal plus noise, and Фф, (t) is the 
crosscorrelation function of the signal plus noise, and the signal. 

The requirements that must be fulfilled to use Eq. 2.1 are 
(1) that all signals are stationary and exist for all time; (2) that 
the signal and the noise are uncorrelated. In order to solve the 
equation, the autocorrelation function of the signal and the auto- 
correlation function of the noise must be known. These may be 
obtained by taking the Fourier transform of the respective power 


ж 
spectral densities, as given by the Wiener-Khintchine equations: 


e 


NO EEE [ne "er om [es | 


- 


"sepe C. Newton, Jr., Leonard A. Gould, and James F. Kaiser, 


Analytical Design of Linear Feedback Controls (New York: John Wiley 
and Some, Inc., 1957) p. 144. 


4 
"John L. Stewart, Fundamentals of Signal Theory (New York: 
McGraw-Hill Book Company, Inc., 1960) p. 306. 
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со 
1 ών ο Doe 
P&T) = OS dw = 2n 4 ЁС 
= co 
where М (Ш) is the power density spectrum of x(t). Even when 
sufficient knowledge of the signal is available, the solution of 


Eq. 2.1 is generally an extremely difficult task, 
1. SOLUTION BY SPECTRUM FACTORIZATION 


The solution of Eq. 2.1 using the technique known as spectrum 


2١.97 The 


factorization has been developed by several authors. 
autocorrelation function of the signal plus noise, ф,,(8), 15 


factored so that 


4 = 
PSF PF E) ¥ E) со 


where CY contains those factors of Po (8) which have poles and 


zeros in the left half plane. The term 


TD (2.1-2a) 


ф MC М 


is defined by those terms of the partial fraction expansion of: 


Q (s) 


ZX 


Q (8) 


ZZ 


which have poles in the left half of the s-plane. That is 





φ p (8) q (s) 
--- -|-ᾱς-- + | ته‎ (2.1-2b) 
Ф 22 i ee е M zz 9! 
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Ihe resulting weighting or transfer function of the filter is given 


by: 


(2.1-3) 


The use of this equation may be demonstrated by the example 


shown in Figure 1. Consider that x(t) has an autocorrelation 


function: 
o (T) = al ад (2.1-4) 
XX 2 қ 
Тһе Laplace transform of the autocorrelation function is then 
φ. (5) - — (2.1-5) 
ΧΧ 2 ١ 
1 - 8 
which yields the power density spectrum, 
1 1 
К (w) = = — (2.1-6) 
x 217 1 + 42 


This signal can be shown to correspond to the random 
walk signal shown in Figure 2, where x(t) is a rectangular wave with 
values +; ΟΥ - 1 and with zero crossings located at event points 
which are Poisson distributed with an average frequency of 5 crossings, 


per second, This time domain signal is not unique since many signals 


may have the same power density spectrum. 


Now consider a white noise source where, 
факт). 2 507) (2.1-7) 


where $ (т) is the Dirac delta function, which has the following Laplace 
transform: 


Ф 65) = 1 (2.1-8) 


l6 


TIME —————J 


FIGURE. 


RANDOM WALK SIGNAL 
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The autocorrelation function of the signal plus noise is 
= F + ar = 
o) re 212 2 с (2.1-9) 


Since the signal and noise are assumed to be uncorrelated the last 


two terms are zero and q, , (S) reduces to: 


psp (в) 
2 
1 —— ТЕ o (2.1-10) 
1-5 1 - 8 


The crosscorrelation function of the signal plus noise, and the signal 


is: 
PE) = / (v(t) + x(t)) x(t 3 T)dt 
= | x(t +T)dt + ES πι ШЕ 
= | x(t) x(t * T)dt - P&T) (2. IE 


- 09 


Since the noise and the signal are uncorrelated, 


= зо дав 
Px 9) 7 2 
1 - s 


(21-120 


Factoring Eq. 2.1-10 results in: 


— 


2 
9 Gy. ale o a алан 


1 - Ж; (1 + s) s» 


Therefore, from Eq. 2.1-1, 
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+ 
$ 
+ 
0 


Ф ير‎ A 
and, 
πο " i - Ξ 
Νου, 
9 ааай 1 
Ф „„(в) аас =) ως αμα, 
1 1 1 


Р (2.1- 14) 


2) 14 ME S 1 + 5 


Therefore, from Eq. 2.1-2a, 


1 1 (5) يم 
a = 2. ia E | + " (2.1-15)‏ 
φ ЖЕ.‏ 

+ 


Substituting into Eq. 2.1-3, 








H(s) à s 
2.208 ILF s V2 + s 


l l 
2.41% | +в (2.1-16) 





The resulting filter weighting function, or impulse response, is: 


1 - 2 є a. 


wa E > 





H(t) = 
0 ut) (2.1-17) 


This filter is physically realizable, since it is causal. 
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2. THE BODE-SHANNON SOLUTION 


The frequency domain technique for the Wiener-Kolmogorov filter 
as developed by Bode and dT is to convert the signal power 
spectral density to that of white noise, and then have a filter 
operate on this white noise, The development of this technique 


А 
του and is discussed in Appendix 


appears in several references 
B. The following results are obtained, 

Given WW), the power spectral qes s of the signal x(t); 
and πώ), the power spectral density of the noise v(t), consider 


the transfer function of the following minimum phase (zeros and 


poles in the left half of the s-plane) filter: 


Bw) = VW @) +W (0) (2.2 


From a white noise input this filter will produce a power spectral 
density equal to that of the signal plus noise. As shown in Appendix 
B, the following transfer function can then be obtained for the non- 
causal optimal filter. 


WW) B(w) 


тр - TW O) FW DJ )2.2-2( 


The preceding equations are represented in Figure 3. Since h, @) is 


non-causal, let 


Ralph .ل‎ Schwarz and Bernard Friedland, Linear Systems 
(New York: McGraw-Hill Book Company, 1965), p. 334-45. 
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FIGURE 3 


BODE-SHANNON THEORY 
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0 Ἵππον. 
g, (t) - (2.2-3) 
> 
h (t) Сг 0 


where g(t) represents the portion of h (t) which exists for t 2 0 
This is equivalent to the spectral factorization of Eq. 2.1-1 and 
Еа. 2.1-2а. Тһе desired optimum causal filter has the transfer 


function 


H(w) = — (2.2-4) 


Now consider the example problem again. The signal, x(t), has 


a power spectral density 


1 1 
М (0) 5 = — (2,2-5) 
x ОТЕ 1 + 02 


апа the uncorrelated measurement noise is white, with a power spectral 


density 


mod 3 
Ww) zi (2.2-6) 


The power spectral density of the signal plus noise is therefore 


2 
WE) н у зи уа —— πι. 22-8. (2.2 
j : х 2π |l + 2n |l +O 


This may be factored 


BR Va + 19) (V2 - jw) ч 
о) ко (ἜΤΙ, (1 - jw) (2.2-8) 


Thus, the minimum phase transfer function required to produce this 


spectrum from white noise, using Eq. 2.2-1, is: 
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(2.2-9) 


The inverse of this transfer function will produce white noise from 
the signal plus noise. 


Now, substituting in Eq. 2.2-2, 


I 


һ Сш) 





N E or 
(2 + ш2у (14 jo Van 


1 1 


(V2 - jw) (1 + jw) Von 


κ“ „= :.. τω (2.2-10) 
2.414 Y 2π o» jw 1 + jw 


Е 
= 


By taking the inverse Fourier transform of Eq. 2.2-10, the impulse 


response, or weighting function, is then 








ҰЛТТЫ د‎ E 
b, C? - 
1 7 
F š = ЕТЕ тг )2.2-11( 


which is non-causal. The best causal impulse response, from Eq. 2.2-3, 


15 
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gj (t) 7 
1 1 -t 
2 
Моп 2.414 € анг DB) 
The desired impulse response is therefore, 
G (w) 
ü) 
H(w) = 8 (9) 
/ \ d 
= >” |, ен. | 1 + ju + ju V 2n 
IT 12.414 Ut», 4? + jo 
le 1 1 
2.414 J2 * ys 
and 
0 12-00 
H(t) = 
: „u i tz0 (2.2-13) 


2.414 ; 


This filter impulse response is identical to Eq. 2.1-17, as expected. 


3. DISCRETE TIME SOLUTION OF THE 


WIENER-KOLMOGOROV FILTER 


Ihe initial problem to be solved for digital computer simulation 
of the Wiener-Kolmogorov filter is the realization of the discrete 


signal, x(kT), based on the required power spectral density 





(2.3-1) 


жау = 2 


1. 

2π Іш 
Since the power density spectrum of white noise is a constant, this 
power spectral density may be obtained by driving a linear system, 


G (jo), by white noise, w(t), so that 
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πω) =G (w) GC) w у= -- |1 (2.3-2) 
211 1 + ىن‎ 


where 


i E 1 
G, (30) 2221 4 po 


It is extremely difficult, if not impossible, to obtain a 
discrete white noise source which is flat over all frequency. Іп 
practice, what needs to be done is to make the spectrum of the noise 
source flat over the bandwidth of G (jo). The 3 decibel cutoff 
frequency for G, (30) is one radian per second. The effective noise 


bandpass is given by 


B -5 (1) = radian/sec. = 1/4 Hz. (2.3-3) 


where noise bandwidth is 2 times the 3 decibel cutoff — 
The sampling rate for the discrete noise should be much greater than 
1/4 Hz., and 100 Hz. was chosen as convenient. The sampling period, 
T, therefore is 0.01 seconds. The spectrum of the noise should then 
be flat to 50 Hz. as shown in Figure 4. 


The autocorrelation function of this band limited white noise 


is then given by 


100 sin (100 TT 


Py T2 100 πτ Vo MD 


as shown in Figure 5. The variance of w(t), the noise sequence, is 


given by 


(2.3-5) 


Э [= 


= 2 _ = 
ο σα 100 


"Mischa Schwartz, Information Transmission, Modulation, and 


Noise (New York: McGraw-Hill Book Company, Inc., 1959), pp. 206-207. 
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FIGURE 5 


AUTOCORRELATION OF WHITE NOISE SOURCE 
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Ihe white noise source is then simulated by a series of pulses of 
repetition period .01 seconds, whose аш со have a Gaussian 
distribution with zero mean and variance 100, followed by a zero 
order hold as shown in Figure 6. The recursive equation necessary 
for the computer simulation of the desired signal source is derived 


by taking the Z transform of the input-output transfer function of 





Figure 6. 
/ 
x(2) _, 1 і d 
w(z) s + 1 s | 
| 
-Т 
s — (2.350) 
27-06 

Непсе 

х(2) (2 - ee = w(z) (1 - D 
and 


zx(z) = x E + w(z) (1 - = 


Therefore, 


A =e) =) aie aan oe 0: 


To check that the signal does, in fact, have the desired power 
Spectral density, a subroutine entitled HARM” was used to obtain the 
Fast Fourier Transform. The Fourier coefficients were then squared 
and plotted versus radian frequency, and the graph, Figure 7, obtained. 


Included on this graph is the desired power spectral density curve. 


μα. πη Business Machines Corporation, System /360 
Scientific Subroutine Package, (360A-CM-03X) Version II Programmer's 


Manual, (White Plains, New York: International Business Machines 
Corporation, 1967) 
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G16 


n 
512 


316 


The program was run for 2048 samples and it can be seen that the 
power density spectrum approximates the desired spectrum. The 
computer program for generating the signal and checking its power 
density spectrum is included as Appendix C. 

Now that the proper signal has been obtained, a recursive 
equation for the filter is needed. The filtering system is drawn 
in Figure 8. Using the Z transform, the recursive filter equation 


is developed as follows: 








Hence | 
3(2) E - SVP 2 - TAE 1 - πα z(z) (2.3-9) 
and 


zQ(z) = (2) e م‎ Z(k) 


о о 


V2 T 1 ΠΠ 


(23-10) 


Therefore, the recursive equation is: 


Z(k) (2.3-11) 


So + 1) - wa Ў (к) + 2417 E - e V2 T 





The error is 


E(k) = Хас) - X(k) (2,351 
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ЖИЕ 


НА1114 AOHODONIOM-UHNSIA 2225 


3 م 





A 
The assumption is made that X(0) is zero. The sampling interval 
is 0.01 seconds. The measurement noise, v(t), is white, with a power 


=, and is developed in an identical fashion 


spectral density, 1212 Ξ 
to w(t), the driving noise for the signal generator. The two noise 


Sources are uncorrelated. The mean square error and the variance 


of the error are then calculated, on a continuing basis, using the 











equations: 
1 n 
E(n) 2 X  E(k) (2.3215) 
kel 
n 
ο τα πο. (2.3310) 
п kel 
2 
с = ΙΙ, EO (2.3-15) 


Ihe computer run was for 30 seconds, i.e., 3000 samples, so that 
steady state is reached. The results are discussed in Chapter V. 
The computer program for the Wiener-Kolmogorov filter appears as 


Appendix D. 
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CHAPTER III 
BODE-SHANNON DISCRETE TIME FILTERING 


The Bode and Shannon frequency domain approach may be applied 
in the time domain directly, as developed by Schwarz and BELAM a l 
The formulation is now in terms of summations, discrete signals, and 
discrete power spectral densities. This technique results in the 


following equation which is similar to the Wiener-Hopf integral 


equation. 


Pay (kon) = h(n, i) p, (j,k) (3. 1a) 


ΤΊ la 


=0 
k = 079 эт δὲ... m" 
where 


P (k,n) =E E > (3.15) 


is the discrete crosscorrelation function between z(k) and x(n), 


and 


p. k) = E m 22 08-16) 


is the discrete autocorrelation function of z(k). The unit impulse 
response of the filter is h(n,j). The notation h(n,j) denotes the 
response at time n due to an impulse applied at time j where j = n. 
Since the impulse response of the filter must be causal, this requires 
that h(n,j) equal zero for n less than j. Equation 3.1 holds. for both 
stationary and non-stationary signals. When the signals are stationary 
then 


p,„(Kın) = p. n - k) (3624) 


P (k) = PŒ- j) (3.2b) 
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Apta = J) (3.26) 
and Eq. 3.1 reduces to 


q,,) 7 E Πα) 9, (m- i) (8824) 


m=0 


where 


P. 
il 


n-k andm 7n - j (3.2e) 


The power spectral density of a discrete signal, z(kT), is defined 


х 
as: 


W_(z) 


TZ [Ж 


ΤΣ φωτ) α΄" (3.3) 


1 =-0 


This is the discrete form of one of the Wiener-Khintchine equations. 
The prime denotes the power spectral density of a discrete signal. 

It should be noted that the z variable within the parenthesis denotes 
the Z-transform shifting variable. Proceeding as in Appendix B, a 
transfer function, B(z), is developed which will convert a signal with 
a power spectral density of W, (z) to a white noise process, u(kT). 


This requires that 
m ! 
B(z) B(2 ) W (z) = 1 (Ύ 


Since B(z) is causal, its Z transform must be given by 
co 


B(z) = У b(n) z^ (3259 
п=0 


4 
"Schwarz and Friedland!” use this definition. Kuo uses a factor 
ОГТ instead of T. 
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1 
τε W (2) can be factored into a product, 





W (2) - W (z) W (271) (3.6) 
then 
Е (3.7) 
0 


For a white noise input ΕΩΣ equals 6(j,k), and Eq. 3.1 reduces to 
Px (kom) = h@, j) ΤΙ 
p Z (on) = 0 for j 2 n (3.8) 


The impulse response of an optimal causal filter for a white noise 


input, u(kT), and a desired output, X(kT), is then given by 
g(n,j) = 9 (k,n) (3.9) 


For stationary processes Eq. 3.9 reduces to 


в(т) = Ф (в) for m = m - k > O (3.10) 


ux 
Taking the Z transform of Eq. 3.10, with the aid of Eq. 3.3, yields 


| w (z) 
g(z) =- EE (3.11) 


C 


١ 
where W | is the cross power spectral density of the white noise, u(kT), 


and the signal. The subscript c denotes the causal part of the transform 


Р -1 
(terms involving z ). However 


ΩΩ, - 7 (а) в(27 +) (355 5) 


so that 
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| 5 
МЕС) 
g(z) = — 2222 (ӨГІЗ) 


The optimum filter is thus 


H(z) 


G(z) B(z) 


и, (2) CR 


= B(z) (3.14) 


H 


The procedure of taking "the causal part of" is analogous to the 
EN 


Spectrum factorization done in Chapter II. 

Using the above technique, the previous example is now solved 
for discrete time signals. The signal and noise correlation functions 
are 


1 -|k|T 
Page (KT) 2 Ё 


ll 


(3.15) 


ll 


Ф (KT) = 6 (k) (3.16) 


The sampled power spectral density of the signal may be obtained by 


means of Z transforms and is as follows: 


Т 2 3 In 


1 - << 
= T n е (3.12 


(l= eg) (1-0) 


! 


W, (x) 


where € equals an For the noise, the sampled power spectral density is 


١ 
W (z) = T (3.18) 
ν 
Therefore, the power spectral density of the signal plus noise input 


to the filter is 


ll 


W (z) W (2) +W (z) 


a (3.19) 


я Т 21 
)1 - €z ( )1 = ez) 
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This may be written in the following manner for ease in factoring: 


1 


ej- 2E 
Ч) b (1 - ez 1) (1 - ez) (3.20) 
where 

b s Ā — - 0031221 


The transfer function needed to produce white noise from the signal 


and noise input to the filter is then, 





| -1 
B(z) = Мъ/ет de |- ل‎ (3.22) 
1 - 2 W, (z) 
Since the signal and the noise are uncorrelated, 
p Ü (k) = Ф, (К) (3.23) 
and 
! ! : 
ч) 7 W (2) (3.24) 
ius, from Eq. 3.134 
| -1 
W (z) B(z ) 
(2) = шэг. 
8 T 
ο 
1 - c? Cu zl ШЖА- 62) 
= FEST OO A b/eT . e a 03:25) 
(1 - ez ) )1 = 62( (1 - bz) i 
c 
Simplifying and taking the causal part yields 
2 : 
1 - © 1 
EZ) μμ η а] (2526) 


The discrete time transfer function for the optimum filter is then 


H(z) = g(z) B(z) 
a DE | c U . (3727) 
Te 11 - δε 1 Ж ἆ 
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Ihe optimum filter diagram is shown in Figure 9. Since b, €, 


and T are constants, H(z) may be rewritten as: 


H(z) =C : 1 (3.28) 
1 - bz 
where 
2 
Ъ(1 - e^) 
Ç Te(l - be) (3.29) 
The recursive equation for the filter is developed as follows: 
Х(2) C 
R(z).- ον © 
1 - ра 
Z (z) (1 - bz >) = Cz (z) 
-1 A 
%(z) = bz ` x(z) + Cz(z) 
#(к + 1) = DÉ(K) + CZ(k + 1) (3.30) 
The error is 
4 
E(k) = X(k) - X(k) (3.31) 


For the specific example the coefficients for the filter, Eq. 3.30, 
become b = .869 and c = 12.45. These differ from the coefficients for 
the discrete Wiener-Kolmogorov filter, Eq. 2.3-11. One might expect 
that these results would be identical, but although Schwarz and ME S ' 


consider both developments the discrepancy is not explained and remains 


a subject for further investigation. 
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CHAPTER IV 
THE KALMAN FILTER 
1. DISCRETE TIME SOLUTION 


The objective in this chapter is to examine the filtering 
problem from the state point of view as developed by inan ome 
In the previous cases, the signal was characterized by either the 
autocorrelation function or the power spectral density. Here the 
signal is characterized by the output of a linear dynamic system 
driven by white noise similar to Bode-Shannon. This presents no 
problem, for if the power spectral density is rational, it can 
always be represented as the output of a linear system driven by 


white noise. As in Chapter 11, the linear system has the transfer 


function, 


Ж ρε 
G) (s) = sasa) (4.1-1) 





Since this is a first order system, the observation matrix is unity 
and scalar notation can be used. The corresponding differential 
equation is 


ante) = E) w(t) (4.1-2) 


Then, from state variable theory, 


Š = zd and ne) E 


and the differential equation may be expressed in equivalent recursive 


form: 


ока) = е у) (4.155) 
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The signal modelis outlined in diagram form in Figure 10; w(t) 
and v(t) are independent gaussian random processes with zero means 


and covariances given by: 


cov | «(ЕЭУуеч (139 = M8 (£ - 1) 
cov [v(t), v(T)] =RÓ (Е - т) 
cov [w(t), v(T)] = 0 


Q and R are the variances of the signal white noise source and the 
measurement noise, which are set equal to unity in order to duplicate 
the model for the Wiener filter. The noise sources are simulated in 
the exact manner as described in Chapter 11, Section 4. 
| | . А | : : 
the Kalman filter is shown in Figure 11. G(k) is an adjustable 


gain which minimizes the mean square error, 





BE (2) = 


5 ir 
IR 
| πα = 


(E(k))? (4,1-4) 
1 


where 
E(k) = X(k) - X(k) 


G(k) is determined by the following equations: 


PUK) - (1 - οἷο POETE 21 
πι Πως ρω hp CE 


where P(k/k - 1) denotes P at time k given the value at time k - l. 


4] 
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ЖЕЗ 


T3G0N 1949215 غ117‎ 


ОТ 3 


(т +з)х ο. - (3)A 


Kered 


(A)X 


(3)4 
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NVWIVM HLHdHdOSIG 


11 3 


21343113 





1 2 
[2(а+1) - 2(а+1/п) | 
0 


{к 
o 
| 


в(к/к-і) “в 0 (209-2(к/к-1)7 


! M 


11 


2 


II 
Im 
ЕТ ас 


P(k/k) 5 E [zao - IJ £ E - 22 


нэ, 
! 
سم‎ 


also, 


It is necessary to define P(1/0) in order to start the recursive 
process. 

It is important to note that once P(1/0) is specified the gain 
equations, Eqs. 4.1-5, 4.1-6, and 4.1-7, are recursive and do not 
depend on the observations Z(kT). They depend upon the variance 
of the signal driving noise, the variance of the measurement noise, 
and the characteristics of the linear system used for signal generation. 

From Figure 1], the recursive equation for the Kalman filter may 


be expressed as 
X(k) - $X(k - 1) * G(k) | σα - 8X(k - D | (4. 1-8) 


The error is 


E(k) = X(k) - X(W (4,1-9) 


As in the Wiener simulation, the sampling interval is 0.01 
seconds, $ (0) is the initial estimate at time t = O and the length 
of the computer run is 3,000 samples. The computer program for 


the Kalman filter simulation is included as Appendix E. 
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2. EQUIVALENCE OF KALMAN AND WIENER-KOLMOGOROV FILTERS 


IN THE STEADY STATE 


6 | 
It has been shown by C, E. Hutchinson that the continuous 
Kalman and Wiener filters are equivalent in the scalar case after 
the steady state condition is reached, i.e., G(k) is a constant. 


This equivalence is developed as follows. 
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Given a message spectrum, WW) = 22 
bo +1 


2 
a noise spectrum, М = c 


and a crosscorrelation function, © (1) = 0 


The optimal causal filter, in the Wiener sense, is then 


H(s) = (k, - k) -+ (4.2-1) 


2 


where 


1/0 


1/bc V а^ + 22 


^ 
lI 


^ 
u 


Now consider the Kalman filter. The signal model is, in general, 
defined by: 
παω ο) ος) (4,2-2) 


z(t) 


Mx(t) + v(t) (4.2-3) 


The best estimate of x, designated X, is obtained from the continuous 
filter equations which correspond to the discrete equations 4.1-5, 


4.1-6, 4.1-7, and 4.1-8. 
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%(t) = FR(E) + g(t) [α(ε) - Mx(t)] (4.2-4) 
where 

g(t) = P(t) MR 

P(t) = FP(t) + عم‎ - P(t)M R7 MP(t) +Q 
and 

cov [w(t), w(T)] = Q6(t - τὸ 

cov [v(t), у(т)] = R6(t - 7) 

cov [w(t), vr) ama 


If a one dimensional problem is considered with 


Е =-1/Ъ 

Б "me 
R = эж 
M= 1 


the Kalman solution becomes 


II 


Сс) а «1/5 #(t) + g(t) lat) xtt)] (4,2-5) 


where 
g(t) = P(t)/c* 


P(t) du 


ορ ο προ aT 


In the steady state, when P(t) is identically zero, the equation for the 
Kalman filter reduces to 


YO = -1/b x(t) + 


Гал IDE icu qua о ME) | (&.2-6) 
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А 
The transfer function relating the estimated signal, X(s), to 


the signal plus the noise, Z(s), can be expressed as 
X 1 
هی‎ μα k) — (52227) 


which is identical to Eq. 4.2-1. This result can be used as a check 
on computer programs for the Kalman and Wiener-Kolmogorov filters. 
When G(k) has reached a constant value, the mean square errors should 


be identical. 
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CHAPTER V 
RESULTS AND COMPARISONS 


The continuous Kalman and Wiener-Kolmogorov filters have been 
shown, in Chapter IV, to be equivalent in the steady state condition, 
that is, when the Kalman gain is constant. This result is verified 
experimentally, as shown by the computer curves of Figures 16, 20, 24, 
and 28. The steady state gain of the Kalman filter reaches a value 
of .0041335 as expected. Аз soon as the gain reaches its steady 
value, the curves of mean square error versus time (Figures 18, 19, 
22, 23 20, DO and 32) are identical to each other and to the 
Wiener-Kolmogorov filter curves (Figures 14 and 15). 

In order to start the recursive process, the Kalman filter 
requires a priori knowledge of the initial error covariance, P1/0* 
Four different values of 51/0 (0.0, .414065, 1.0, and 10.0) were 
chosen to compare the effectiveness of the Kalman filter for different 
values of thís term. The value of P1/0 = ,415065 was chosen because 
it gives the steady state gain (.0041335) at the start of the computer 


run, Withithisap the performance of the Kalman and Wiener-Kolmogorov 


1/0 
filters are Identical (Figures 13-15 and 21-23). 

In simulating the non-stationary signal, two different initial 
values were chosen; one to correspond to a small initial error in 
estimation and the other to correspond to a large error in the initial 
estimation. In a tracking problem, the first case would represent 
the situation where a target is picked up at long range. The second 


case represents the situation where a target suddenly appears at 


relatively short range. For small initial errors, the Wiener-Kolmogorov 


48 


and the Kalman filters performed essentially the same (Figures 13, 
14, 17, 18, 21, 22, 25, 26, 9, апа 30). ог large initial errors 
the Kalman filter, with large initial covariance of error, Р уо» 
gives the best transient response (Fig. 31). For large initial 


error, if P is equal to zero, the transient response of the 


1/0 
Kalman filter takes somewhat longer time to settle than the Wiener- 


Kolmogorov filter as shown in Fig. 19. 
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CHAPTER VI 
CONCLUSIONS 


The responses of the discrete Kalman filter for various initial 
error covariance and the discrete Wiener-Kolmogorov filter have been 
compared for a non-stationary random input with defined statistical 
properties. It was necessary to develop a computer program for the 
simulation of the desired spectrum from a white noise source. The 
discrete time Wiener-Kolmogorov filter was also programmed for the 
digital computer. 

It has been shown that when the gain of the Kalman filter is 
equal to its steady state value the performance of the Wiener- 
Kolmogorov and Kalman filters are identical. For large initial error 
the Kalman filter (with a large initial gain) shows the best transient 
response. When the initial error is small the Wiener-Kolmogorov and 


the Kalman filters respond in as essentially equivalent fashion. 
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APPENDIX A 
DEVELOPMENT OF THE WIENER-HOPF INTEGRAL EQUATION 


Ihe Wiener-Hopf equation gives the relationship between the 
filter weighting function, the autocorrelation function of the 
input to the filter, and the crosscorelation function of the filter 
input with the desired filter output. The problem is described in 

š : : š š ; A ; 
Figure 1; x(t) is some signal, v(t) is white noise, x(t) is the 
estimate of the signal, and e(t) is the error in estimation; reduce 
e(t) to a minimum in the mean squared sense by a correct choice of 


h(t), the impulse response of the filter. The error by definition is: 


e(t) 2 3D - &(t) (A.1) 


Squaring both sides yields, 
ο = ас - 22(t) x(t) 1-7 (A.2) 


The output, x(t), and the impulse response of the filter, H(t), may 


be related by means of the convolution integral. 


со 


R(t) = Ju x(t - T) dr (A.3) 
=% 
where T is the variable of integration. Squaring Eq. A.3 yields 
со 
τα = Jun z(t - T)dr Γκαρ z(t - EP (A.4) 
= GO _ 


Substituting Eq. A.4 into Eq. A.2 gives, 


e2(t) = x(t) - 2 2 


Pato) z(t - T)dr HQ) z(t - TdT, (A.5) 


NC 
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By definition, the mean square error is, 


1 


1 
lim 27 / e? (t)dt (A.6) 
Tro “τ 





e^ (t) 


Substituting Eq. A.5 into Eq. A.6 yields: 


со T со 
2 ) 1 2 Р 1 
e (t) = lim η x (t)dt - 2lim oT dt H(T)m(Ct == (Ct dt 

τὸ EN Too _T ls 

1 со со 

zu 
+ Lim ir | 4 June - T)dt | H(t, )2¢t - Tat, (A.7) 
-1 =o - 6 


The definition of the input autocorrelation function is 


1 


9, , CT) = Lim > / z(t)z(t ^ T)dt (A.8) 
xm 


Similarly, the crosscorrelation function between the input and the 


output is 
T 


9, 4T) 2 „Lim 5, | RCE + Det (A.9) 
Е 


and the autocorrelation function of the output is, 
^ T 
1 A A 
Т = 1 То zi А, 
92400 = dim d, | ote Tat (A. 10) 
-T 
Now, making use of the correlation functions, and interchanging the 


order of integration in Eq. A.7 yields, 





= = Цэ 09) - 2 J H(T) e, , C dT ES 
| нег [rev та 017 (А.11) 
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To determine the system function that minimizes the mean square 
error, assume that a solution exists and denote this solution as 


H (t). Now construct a weighting function: 


H(t) = H (t) + ôH, (t) (А.12) 


where H (t) is the assumed solution, H, (t) is any arbitrary realizable 
weighting function, and ὃ is a parameter which may be varied to test 
whether H (t) is the solution. The construction is such that at Š 
equal to zero, the derivative of mean square error with respect to Š 
must be equal to zero. By setting this derivative equal to zero for 
6 equal to zero, a condition for H (t) which must be satisfied in 
order for it to be a solution,is specified. 
Substituting H(t) as given by Eq. A.12 into Eq. А.11 and differenti- 


ating with respect to Ó yields: 


d Б Ë 


45 = - 2 | Hg (TI, (TdT + | & (ar | H.G 9, „(Т - TdT, 





4 EX ar | H (T) q a Cr - өн) 
* ЯГ (T)dT 1 H. (T) фат - TT, (A.13) 


Because of the even property of autocorrelation functions of stationary 


signals, 


р (Т,-1) =p_ (T- T) (A.14) 


zz 1 22 1 


which leads to 


25 


с со 


une | Hs (Ti) pan r 3 TT, 


= [gar аа πια, (Α. 15) 


Substituting Eq. A.15 in Eq. A.13 and setting 





Ξ 0 at 6 20 


dó 
yields: 
с с 
2 | EG E OO U CO 
~ СО = 


Since Hs (T) is a realizable weighting function, it must be zero for 
values of t less than zero. The only way that Eq. A.16 can be 
satisfied for Т 2 0 is that the factor in the brackets be equal to 


zero. Thus, 


eo 


Їх (71) NT - TdT, = DN for T = 0 (А.17) 
-00 
Equation A.17 is the famous Wiener-Hopf integral equation. 

Strictly speaking, the H (7) that satisfies Eq. A.17 produces a 
stationary value of the mean square error and does not necessarily 
ensure a minimum. It can be shown that the solution of Eq. A.17 does 
yield a minimum, by considering the second derivative of the mean 
square error with respect to 6. Differentiating both sides of Eq. A.13 
yields 


= 3! H, (T)dT | Hs (T4) ШЕ - TdT, 


B 
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This is the same as twice the mean square value of the noisy input 
signal filtered by a weighting function, Н, (1). Since this mean 
square signal can never be negative, this shows that the second 
derivative of the mean square error is always positive, hence the 


solution of Eq. A.17 corresponds to a minimum. 
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APPENDIX B 


DEVELOPMENT OF THE BODE-SHANNON SOLUTION 


Bode and Shannon used frequency domain techniques to solve the 
filter problem by converting the actual signal power spectral density 
to that of white noise, and then operating on this white noise signal. 
The solution may be described by considering Figure 1; x(t) is the 
signal, v(t) is white noise, $(t) is the estimate of the signal, and 
e(t) is the error. The assumptions are that, (1) input signal and 
noise are uncorrelated and stationary, and (2) all time functions are 
Fourier transformable. 

The impulse response of the filter, h(t), is chosen on the basis 


of minimum mean square error. The error signal is defined as 


e(t) = (t) - x(t) (B.1) 


and its Fourier transform is, 


E(w) Ж (ш) - X(w) 


(ха) n νω)] h(w) - X(w) 


Гао) - 1] X()  πωγνώ) (B.2) 


Since the signal and noise are uncorrelated, the power spectral density 


of the error is 


2 
w (w) = ће син i м. (ш) + | (9) | Ж (B.3) 


Hence the mean square error is 
со 


e^ (t) = — J ES - 11% 1,60) + [hw] w w) 


- 


do  (B.4) 


— 


This is to be minimized by the proper choice of h(%). 
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Note that to obtain a solution, the power spectral densities 


ο and κ) must be specified. Let 


C (9e 399) 


h(w) 


C(w) (Сов0(0) + j sin9(w)) (B.5) 


Substituting into Eq. B.1 and simplifying yields, 
e 


ee) => КО; + ζω + 1 - 2C(0)cos(9(0)) P u, Ww) du 


-οὉ 





(8. 6) 


Since C(w), W. (Ὁ), and W (0) are nonnegative this expression is 
minimized when cos 9(0) has its largest value; that is, when 09(w) 


equals zero. Thus, 


2 "uu 2 
e apes") = = | E (9) 





нө» + ὦ] - 2C(9)W. (9) 4221 аш 





(В.7) 


In order to find the CW) that minimizes this expression, complete 


the square of the terms in the bracket. Thus, 





2 = 2s 
E oa? ~ 2m 
W (w) 7 ы (шн (0) 
c(w) VW (w) + м, (ш - E ge n 
E 
v WW) ἘΝ Ὁ) W. @) ни ὦ) 


(В.8) 
Eq. B.8 is minimized when 
W (w) 
ο эвт (8.9) 


W. (9) ἘΝ Ὁ) 
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with minimum mean square error given by 





W (@)W (w) 
2 Р! x V 
كن‎ СЕ) E or — rn d (B.10) 


W. (2) 0 (9) 


Since 0(@) equals zero, 
ШКС» 


h(w) = C(w) = --------- (B.11) 
Ν (9) رم واي‎ 


The impulse response is therefore given by 


м (о) 


іше 
W (w) +W (w) 
x V 


h(t) = e dw (B.12) 


2TT 


ζω 
which, in general, does not vanish for t « O. Thus Eq. B.11 is 
generally not physically realizable. 

Ihe minimum phase filter which converts a white noise input to a 
spectral density equal to that of the signal plus noise, must have 


a transfer function of magnitude 
- w) +W (wW 
B (w) vn ) 2 ) 


Since Eq. В.11 gives the best noncausal operation on the signal plus 


noise, the best noncausal operation on white noise inputs requires that, 


W. (@) V Wœ) ا نا‎ 


H, (W) = ey) FW) (B.13) 


This however, is still noncausal. The causal filter is then constructed 


as follows: 


g (t) = 


h, (t) М! (B.14) 


80 


Ihe desired causal filter for operation on signal plus noise 


inputs then has a transfer function given by: 


G (ш) 1 


но) ЕЕЕ ا‎ ------------. 
шин гс) 2 Ум. (0) + WW) 


(B.15) 
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